This notebook contains the code used to generate the figures for the above paper.
In [1]:
from __future__ import division
# Standard library
from collections import OrderedDict
import cPickle
import os
import time
# Third-party
import astropy.units as u
from astropy.constants import G
import emcee
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import root, brentq
from scipy.stats import truncnorm, uniform, scoreatpercentile
from scipy.misc import logsumexp
import triangle
# TODO: check in a stylesheet, use that
plt.style.use("matplotlibrc")
%matplotlib inline
In [55]:
# Number of data points to generate
ndata = 100
# Neutron star distribution properties (fixed)
bounds_NS = (1.3, 2.) # Msun
mean_NS = 1.4 # Msun
stddev_NS = 0.05 # Msun
# Number of steps to use in numerical integration below
Nintegrate = 4096
We generate data with three different models for the true distribution of companion masses:
In [3]:
def mass_function(m1, m2, sini):
return (m2*sini)**3 / (m1 + m2)**2
def make_datasets(N, seed=None):
np.random.seed(seed)
# container to store the various fake data
datasets = OrderedDict()
# container used to store parameters used to generate the datasets
params = OrderedDict()
# ------------------------------------------------------------------------------------------
# Test set 1: a single, truncated Gaussian
#
# generate random primary masses and inclination angles
m1 = np.random.uniform(0.2, 0.4, size=N) # Msun
sini = np.sin(np.arccos(1 - np.random.uniform(0.,1.,size=N))) # uniform in cos(i)
p = dict()
d = dict()
p['name'] = "Gaussian WD"
p['bounds_WD'] = (0.2, 1.44) # Msun
p['mean_WD'] = 0.7 # Msun
p['stddev_WD'] = 0.2 # Msun
p['f_NS'] = 0.
a = (p['bounds_WD'][0] - p['mean_WD']) / p['stddev_WD']
b = (p['bounds_WD'][1] - p['mean_WD']) / p['stddev_WD']
gaussian = truncnorm(a, b, loc=p['mean_WD'], scale=p['stddev_WD'])
true_m2s = gaussian.rvs(N)
d['m1'] = m1
d['sini'] = sini
d['mf'] = mass_function(m1, true_m2s, sini)
d['true_m2'] = true_m2s
d['is_NS'] = np.zeros_like(true_m2s).astype(bool)
params['GWD'] = p
datasets['GWD'] = d
# ------------------------------------------------------------------------------------------
# Test set 2: a truncated Gaussian for white dwarfs + truncated Gaussian for Neutron stars
#
# generate random primary masses and inclination angles
m1 = np.random.uniform(0.2, 0.4, size=N) # Msun
sini = np.sin(np.arccos(1 - np.random.uniform(0.,1.,size=N))) # uniform in cos(i)
p = dict()
d = dict()
p['name'] = "Gaussian WD + Gaussian NS"
p['bounds_WD'] = (0.2, 1.44) # Msun
p['mean_WD'] = 0.7 # Msun
p['stddev_WD'] = 0.2 # Msun
p['f_NS'] = 0.1 # 10% neutron stars
N_NS = int(N*p['f_NS'])
# White dwarf companions
a = (p['bounds_WD'][0] - p['mean_WD']) / p['stddev_WD']
b = (p['bounds_WD'][1] - p['mean_WD']) / p['stddev_WD']
gaussian = truncnorm(a, b, loc=p['mean_WD'], scale=p['stddev_WD'])
true_m2s = gaussian.rvs(N-N_NS)
d['is_NS'] = np.zeros(N-N_NS, dtype=bool)
# Neutron stars companions
a = (bounds_NS[0] - mean_NS) / stddev_NS
b = (bounds_NS[1] - mean_NS) / stddev_NS
gaussian = truncnorm(a, b, loc=mean_NS, scale=stddev_NS)
true_m2s = np.append(true_m2s, gaussian.rvs(N_NS))
d['is_NS'] = np.append(d['is_NS'], np.ones(N_NS, dtype=bool))
d['m1'] = m1
d['sini'] = sini
d['mf'] = mass_function(m1, true_m2s, sini)
d['true_m2'] = true_m2s
params['GWD+GNS'] = p
datasets['GWD+GNS'] = d
# ------------------------------------------------------------------------------------------
# Test set 3: uniform + Gaussian at the Neutron star mass
#
# generate random primary masses and inclination angles
m1 = np.random.uniform(0.2, 0.4, size=N) # Msun
sini = np.sin(np.arccos(1 - np.random.uniform(0.,1.,size=N))) # uniform in cos(i)
p = dict()
d = dict()
p['name'] = "Uniform WD + Gaussian NS"
p['bounds_WD'] = (0.2, 1.1) # Msun
p['f_NS'] = 0.1 # 10% neutron stars
N_NS = int(N*p['f_NS'])
# White dwarf companions
true_m2s = np.random.uniform(p['bounds_WD'][0], p['bounds_WD'][1], size=N-N_NS)
d['is_NS'] = np.zeros(N-N_NS, dtype=bool)
# Neutron stars companions
a = (bounds_NS[0] - mean_NS) / stddev_NS
b = (bounds_NS[1] - mean_NS) / stddev_NS
gaussian = truncnorm(a, b, loc=mean_NS, scale=stddev_NS)
true_m2s = np.append(true_m2s, gaussian.rvs(N_NS))
d['is_NS'] = np.append(d['is_NS'], np.ones(N_NS, dtype=bool))
d['m1'] = m1
d['sini'] = sini
d['mf'] = mass_function(m1, true_m2s, sini)
d['true_m2'] = true_m2s
params['UWD+GNS'] = p
datasets['UWD+GNS'] = d
return params, datasets
In [4]:
test_params,test_data = make_datasets(ndata, seed=4)
In [5]:
fig,axes = plt.subplots(1, 3, figsize=(14,5), sharex=True, sharey=True)
bins = np.linspace(0.,2.,25)
for i,(p,d) in enumerate(zip(test_params.values(), test_data.values())):
axes[i].hist(d['true_m2'], bins=bins, color='#777777', edgecolor='#999999')
axes[i].set_title(p['name'], fontsize=22)
axes[i].set_xlabel("True $M_2$ [M$_\odot$]")
#axes[i].yaxis.set_visible(False)
axes[0].set_xlim(0.,2.)
fig.tight_layout()
In [6]:
fig,axes = plt.subplots(1, 3, figsize=(14,5), sharex=True, sharey=True)
bins = np.linspace(0.,1.4,25)
for i,(p,d) in enumerate(zip(test_params.values(), test_data.values())):
axes[i].hist(d['mf'], bins=bins, color='#777777', edgecolor='#999999')
axes[i].set_title(p['name'], fontsize=22)
axes[i].set_xlabel("$m_f$ [M$_\odot$]")
axes[i].yaxis.set_visible(False)
axes[0].set_xlim(0.,1.5)
fig.tight_layout()
With 'data' in hand, we now start building our model. For all cases, we will try modeling the true $M_2$ distribution as a mixture of a truncated normal distributions. Our parameters will be the mean and variance of the white dwarf distribution, and the fraction of neutron stars.
In [7]:
%load_ext cythonmagic
In [8]:
%%cython
import cython
cimport cython
import numpy as np
cimport numpy as np
cdef extern from "math.h":
double sqrt(double x)
double cbrt(double x)
@cython.cdivision(True)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef inline void _integrand_factor(double[::1] m2, double[::1] mf, double[::1] m1, double[:,::1] integrand,
unsigned int nm2, unsigned int ndata):
""" Cython helper to speed up likelihood calculation """
cdef unsigned int i, j
cdef double mtot, cbrt_mtot4, cbrt_mf
for j in range(ndata):
cbrt_mf = cbrt(mf[j])
for i in range(nm2):
mtot = m1[j] + m2[i]
cbrt_mtot4 = cbrt(mtot*mtot*mtot*mtot)
integrand[i,j] = cbrt_mtot4 / cbrt_mf / m2[i] / sqrt(m2[i]*m2[i] - cbrt_mf*cbrt_mf*cbrt_mtot4) / 3.
@cython.cdivision(True)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cpdef integrand_factor(double[::1] m2, double[::1] mf, double[::1] m1, double[:,::1] integrand):
""" Compute the factor multiplying p(M_2|θ) in the integral of Equation XXX in the paper """
cdef unsigned int nm2, ndata
nm2 = m2.shape[0]
ndata = mf.shape[0]
_integrand_factor(m2, mf, m1, integrand, nm2, ndata)
In [9]:
def py_integrand_factor(m2, mf, m1):
""" Compute the factor multiplying p(M_2|θ) in the integral of Equation XXX in the paper """
mtot = m1 + m2
return mtot**(4/3.) * mf**(-1/3.) / m2 / np.sqrt(m2**2 - (mf*mtot**2)**(2/3.)) / 3.
In [10]:
def likelihood(p, mf, m1, bounds_WD, known_m2):
mean_WD,stddev_WD,f_NS = p
m2s = np.linspace(0., 2., Nintegrate)
dm2 = m2s[1] - m2s[0]
integ_fac = np.empty((m2s.shape[0], mf.shape[0]))
integrand_factor(m2s, mf, m1, integ_fac)
# integ_fac = py_integrand_factor(m2s[:,np.newaxis], mf[np.newaxis], m1[np.newaxis])
# White dwarf companion mixture component
lower, upper = bounds_WD
dist_WD = truncnorm((lower - mean_WD) / stddev_WD, (upper - mean_WD) / stddev_WD, loc=mean_WD, scale=stddev_WD)
# Neutron star companion mixture component
lower, upper = bounds_NS
dist_NS = truncnorm((lower - mean_NS) / stddev_NS, (upper - mean_NS) / stddev_NS, loc=mean_NS, scale=stddev_NS)
p_WD = (1-f_NS) * dist_WD.pdf(m2s[:,np.newaxis])
p_NS = f_NS * dist_NS.pdf(m2s[:,np.newaxis])
# Zero out when evaluating outside of allowed bounds (normally NaN)
integ_fac[np.isnan(integ_fac)] = 0.
p_WD[np.isnan(p_WD)] = 0.
p_NS[np.isnan(p_NS)] = 0.
# we approximate the integral using the trapezoidal rule
integrand_WD = p_WD * integ_fac
integrand_NS = p_NS * integ_fac
p_WD = dm2/2. * (integrand_WD[0] + np.sum(2*integrand_WD[1:-1], axis=0) + integrand_WD[-1])
p_NS = dm2/2. * (integrand_NS[0] + np.sum(2*integrand_NS[1:-1], axis=0) + integrand_NS[-1])
not_nan = ~np.isnan(known_m2)
if np.any(not_nan):
p_WD[not_nan] = (1-f_NS) * dist_WD.pdf(known_m2[not_nan])
p_NS[not_nan] = 0.
return np.vstack((p_WD, p_NS))
def ln_likelihood(p, mf, m1, bounds_WD, known_m2):
return np.log(np.sum(likelihood(p, mf, m1, bounds_WD, known_m2), axis=0))
def ln_prior(p):
mu,sigma,f_NS = p
lp = 0.
# uniform prior on neutron star fraction in the range [0,1]
if f_NS < 0. or f_NS > 1.:
return -np.inf
# uniform prior on position of the Gaussian
lp += -np.log(1. - 0.2)
if mu < 0.2 or mu > 1.:
return -np.inf
# p(σ) = 1/σ (0.02 < σ < 2)
lp += -np.log(sigma)
if sigma < 0.02 or sigma > 2.:
return -np.inf
return lp
def ln_posterior(p, mf, m1, *args):
ll = ln_likelihood(p, mf, m1, *args)
lp = ln_prior(p)
if np.any(np.isnan(ll) | np.isinf(ll)):
return -np.inf
if np.any(np.isnan(lp)):
return -np.inf
if np.all(ll == 0.):
return -np.inf
return lp + ll.sum()
In [11]:
# number of parameters
ndim = 3
def run_inference(mf, m1, bounds_WD, nwalkers=32, nburn=200, nsteps=1000, pool=None, known_m2=None):
if known_m2 is None:
known_m2 = np.zeros(len(mf))*np.nan
# initial positions for walkers
p0 = np.zeros((nwalkers,ndim))
p0[:,0] = np.random.normal(0.8, 0.05, size=nwalkers) # mean
p0[:,1] = np.random.normal(0.6, 0.01, size=nwalkers) # stddev
p0[:,2] = np.random.uniform(0., 0.1, size=nwalkers) # f_NS
args = (mf, m1, bounds_WD, known_m2)
sampler = emcee.EnsembleSampler(nwalkers=nwalkers, dim=ndim, lnpostfn=ln_posterior, args=args, pool=pool)
# burn-in
pos,prob,state = sampler.run_mcmc(p0, N=nburn)
bad_ix = sampler.acceptance_fraction < 0.2
if np.any(bad_ix):
print("Some questionable walkers found...")
print("Resampling their positions based on other walkers")
pos[bad_ix] = np.random.normal(np.median(pos[~bad_ix],axis=0),
np.std(pos[~bad_ix],axis=0),
size=(sum(bad_ix),ndim))
sampler.reset() # throw out burn-in samples
pos,prob,state = sampler.run_mcmc(pos, N=nsteps)
return sampler
In [12]:
# For triangle plots below
labels = [r"$\mu$", r"$\sigma$", r"$f_{\rm NS}$"]
extents = [(0.2,1.4), (0.0,2.), (0.,1.)]
In [14]:
from gary.util import get_pool
# caches the samplers in pickle files
samplers = OrderedDict()
for test_name in test_params.keys():
fn = "{}.pickle".format(test_name)
if not os.path.exists(fn):
pool = get_pool(threads=4)
t1 = time.time()
samplers[test_name] = run_inference(test_data[test_name]['mf'], m1=test_data[test_name]['m1'],
bounds_WD=test_params[test_name]['bounds_WD'],
nburn=500, nsteps=1000, pool=pool)
pool.close()
samplers[test_name].pool = None
with open(fn,'w') as f:
cPickle.dump(samplers[test_name], f)
print("{} seconds".format(time.time()-t1))
else:
with open(fn,'r') as f:
samplers[test_name] = cPickle.load(f)
In [15]:
fig = triangle.corner(samplers['GWD'].flatchain, truths=[test_params['GWD']['mean_WD'], test_params['GWD']['stddev_WD'], 0.],
labels=labels, extents=extents)
In [16]:
for i in range(ndim):
plt.figure()
for chain in samplers['GWD'].chain[...,i]:
plt.plot(chain, marker=None, drawstyle='steps', alpha=0.2, color='k');
plt.ylabel(labels[i])
plt.ylim(*extents[i])
In [17]:
truths = [test_params['GWD+GNS']['mean_WD'], test_params['GWD+GNS']['stddev_WD'], test_params['GWD+GNS']['f_NS']]
fig = triangle.corner(samplers['GWD+GNS'].flatchain, truths=truths,
labels=labels, extents=extents)
In [18]:
for i in range(ndim):
plt.figure()
for chain in samplers['GWD+GNS'].chain[...,i]:
plt.plot(chain, marker=None, drawstyle='steps', alpha=0.2, color='k');
plt.ylabel(labels[i])
plt.ylim(*extents[i])
In [19]:
truths = [np.nan, np.nan, test_params['UWD+GNS']['f_NS']]
fig = triangle.corner(samplers['UWD+GNS'].flatchain, truths=truths,
labels=labels, extents=extents)
In [20]:
for i in range(ndim):
plt.figure()
for chain in samplers['UWD+GNS'].chain[...,i]:
plt.plot(chain, marker=None, drawstyle='steps', alpha=0.2, color='k');
plt.ylabel(labels[i])
plt.ylim(*extents[i])
In [21]:
datafile = os.path.abspath("../data/ELM_WD.dat")
names = ["name", "period", "period_err", "K", "K_err", "mf", "mf_err", "m1"]
elm_WD = np.genfromtxt(datafile, dtype=None, usecols=range(len(names)), names=names)[:-6]
len(elm_WD)
Out[21]:
In [22]:
# eclipsing systems with measured m2
known_m2 = np.ones(len(elm_WD))*np.nan
known_m2[elm_WD['name'] == "J0345+1748"] = 0.72
known_m2[elm_WD['name'] == "J0651+2844"] = 0.50
known_m2[elm_WD['name'] == "J0751-0141"] = 0.97
In [23]:
mf = elm_WD['period']*u.day / (2*np.pi*G) * (elm_WD['K']*u.km/u.s)**3
mf = np.array(mf.to(u.Msun).value)
m1 = np.array(elm_WD['m1'])
In [24]:
fig,ax = plt.subplots(1, 1, figsize=(6,5))
bins = np.linspace(0.,1.4,20)
ax.hist(mf, bins=bins, color='#777777', edgecolor='#999999')
ax.set_title('Data')
ax.set_xlabel("$m_f$ [M$_\odot$]")
ax.set_xlim(0.,1.5)
fig.tight_layout()
In [25]:
if not os.path.exists("sampler_real_data.pickle"):
sampler_real_data = run_inference(np.array(mf), m1=np.array(elm_WD['m1']),
bounds_WD=(0.2,1.44), nburn=500, nsteps=2000,
known_m2=known_m2)
with open("sampler_real_data.pickle","w") as f:
cPickle.dump(sampler_real_data, f)
else:
with open("sampler_real_data.pickle","r") as f:
sampler_real_data = cPickle.load(f)
In [26]:
print sampler_real_data.flatchain[sampler_real_data.flatlnprobability.argmax()]
In [28]:
for i in range(2):
low = scoreatpercentile(sampler_real_data.flatchain[:,i], 50) - scoreatpercentile(sampler_real_data.flatchain[:,i], 25)
high = scoreatpercentile(sampler_real_data.flatchain[:,i], 75) - scoreatpercentile(sampler_real_data.flatchain[:,i], 50)
print("{} -{} +{}".format(scoreatpercentile(sampler_real_data.flatchain[:,i], 50), low, high))
In [29]:
fig = triangle.corner(sampler_real_data.flatchain,
labels=labels, extents=[(0,2),(0,2),(0,1)])
In [30]:
for i in range(ndim):
plt.figure()
for chain in sampler_real_data.chain[...,i]:
plt.plot(chain, marker=None, drawstyle='steps', alpha=0.2, color='k');
plt.ylabel(labels[i])
plt.ylim(*extents[i])
In [31]:
datafile = os.path.abspath("../data/PCEB.dat")
pceb_MD = np.genfromtxt(datafile, dtype=None, names=True, missing='-')
ix = np.isfinite(pceb_MD['P_orb']) & np.isfinite(pceb_MD['K']) & np.isfinite(pceb_MD['M_2'])
pceb_MD = pceb_MD[ix]
len(pceb_MD)
Out[31]:
In [32]:
pceb_MD.dtype.names
Out[32]:
In [33]:
mf = pceb_MD['P_orb']*u.hour / (2*np.pi*G) * (pceb_MD['K']*u.km/u.s)**3
mf = np.array(mf.to(u.Msun).value)
m1 = np.array(pceb_MD['M_2']) # yes, this is 1->2
In [34]:
b = numpy.empty(pceb_MD.shape, dtype=pceb_MD.dtype.descr + [('mf',float),('m1',float)])
for name in pceb_MD.dtype.names:
b[name] = pceb_MD[name]
b['mf'] = mf
b['m1'] = b['M_2']
pceb_MD = b
In [35]:
fig,ax = plt.subplots(1, 1, figsize=(6,5))
bins = np.linspace(0.,1.4,20)
ax.hist(mf, bins=bins, color='#777777', edgecolor='#999999')
ax.set_title('Data')
ax.set_xlabel("$m_f$ [M$_\odot$]")
ax.set_xlim(0.,1.5)
fig.tight_layout()
In [36]:
if not os.path.exists("sampler_real_data_MD.pickle"):
sampler_real_data_pceb = run_inference(np.array(mf), m1=np.array(pceb_MD['M_2']),
bounds_WD=(0.2,1.44), nburn=500, nsteps=2000,
known_m2=known_m2)
with open("sampler_real_data_MD.pickle","w") as f:
cPickle.dump(sampler_real_data_pceb, f)
else:
with open("sampler_real_data_MD.pickle","r") as f:
sampler_real_data_pceb = cPickle.load(f)
In [37]:
fig = triangle.corner(sampler_real_data_pceb.flatchain,
labels=labels, extents=[(0,2),(0,2),(0,1)])
In [38]:
for i in range(ndim):
plt.figure()
for chain in sampler_real_data_pceb.chain[...,i]:
plt.plot(chain, marker=None, drawstyle='steps', alpha=0.2, color='k');
plt.ylabel(labels[i])
plt.ylim(*extents[i])
In [56]:
def model(pars_WD, pars_NS):
# truncated normal for WD's
lower,upper = pars_WD['bounds']
mu,sigma = pars_WD['mean'],pars_WD['stddev']
dist1 = truncnorm((lower - mu) / sigma, (upper - mu) / sigma, loc=mu, scale=sigma)
# truncated normal for NS's
lower,upper = pars_NS['bounds']
mu,sigma = pars_NS['mean'],pars_NS['stddev']
dist2 = truncnorm((lower - mu) / sigma, (upper - mu) / sigma, loc=mu, scale=sigma)
return dist1,dist2
In [57]:
def NS_prob(pars_WD, f_NS, mf, m1, known_m2, bounds_WD=(0.2,1.44)):
p = (pars_WD['mean'], pars_WD['stddev'], f_NS)
P_WD,P_NS = likelihood(p, np.array(mf), np.array(m1), bounds_WD, known_m2)
return P_NS / (P_WD + P_NS)
Make same plot for real data, but no histogram in leftmost plot
In [58]:
def find_confidence_interval(x, pdf, confidence_level):
return pdf[pdf > x].sum() - confidence_level
def sigma_contour(x, y, bins=None, ax=None, levels=None, contourf=False, **contour_kwargs):
"""
Create a contour plot from point density (e.g., posterior samples) and set the
levels to contain the specified fractions of total samples (`levels`).
Parameters
----------
xdata : numpy.ndarray
ydata : numpy.ndarray
bins : iterable, int (optional)
Number of bins along both dimensions, number of bins for each dimension,
or the actual bins (a tuple of arrays).
ax : matplotlib.Axes (optional)
If supplied, plot the contour to this axis. Otherwise, open a new figure
levels : iterable (optional)
Levels to draw the contours at.
contourf : bool (optional)
Use `contourf()` instead of `contour()`.
Other parameters
----------------
contour_kwargs : dict
kwargs to be passed to pyplot.contour()
"""
if ax is None:
fig,ax = plt.subplots(1,1)
else:
fig = ax.figure
H, xedges, yedges = np.histogram2d(x, y, bins=bins, normed=True)
nbins_x = len(xedges) - 1
nbins_y = len(yedges) - 1
x_bin_sizes = (xedges[1:] - xedges[:-1]).reshape((1,nbins_x))
y_bin_sizes = (yedges[1:] - yedges[:-1]).reshape((nbins_y,1))
pdf = (H*(x_bin_sizes*y_bin_sizes))
if levels is None:
levels = [0.683, 0.955]
V = []
for level in levels:
z = brentq(find_confidence_interval, 0., 1., args=(pdf, level))
V.append(z)
X, Y = 0.5*(xedges[1:]+xedges[:-1]), 0.5*(yedges[1:]+yedges[:-1])
Z = pdf.T
if contourf:
c = ax.contourf(X, Y, Z, levels=V, origin="lower", **contour_kwargs)
else:
c = ax.contour(X, Y, Z, levels=V, origin="lower", **contour_kwargs)
return fig
def three_panel(axes, data, sampler, plot_samples=False, mu_lim=(0.2,1.2), sigma_lim=(0.,1.),
pt_alpha=0.5, known_m2=None, panel1_norm=1.):
contour_bins = 27
if known_m2 is None:
known_m2 = np.ones(len(data))*np.nan
map_idx = sampler.flatlnprobability.argmax()
mean_WD,stddev_WD,f_NS = sampler.flatchain[map_idx]
print("MAP params: μ={}, σ={}, f_NS={}".format(mean_WD,stddev_WD,f_NS))
MAP_pars_WD = dict(mean=mean_WD, stddev=stddev_WD, bounds=(0.2,1.44))
MAP_f_NS = f_NS
pars_NS = dict(mean=mean_NS, stddev=stddev_NS, bounds=bounds_NS)
# First panel
m2_grid = np.linspace(0., 2., 256)
if plot_samples:
for i in range(50):
ix = np.random.randint(len(sampler.flatchain))
pars_WD = dict(mean=sampler.flatchain[ix,0], stddev=sampler.flatchain[ix,1], bounds=(0.2,1.44))
f_NS = sampler.flatchain[ix,2]
d1,d2 = model(pars_WD=pars_WD, pars_NS=pars_NS)
func = (1-f_NS)*d1.pdf(m2_grid) + f_NS*d2.pdf(m2_grid)
axes[0].plot(m2_grid, func*panel1_norm, marker=None, color='k', lw=1, alpha=0.1)
d1,d2 = model(pars_WD=MAP_pars_WD, pars_NS=pars_NS)
func = (1-MAP_f_NS)*d1.pdf(m2_grid) + MAP_f_NS*d2.pdf(m2_grid)
axes[0].plot(m2_grid, func*panel1_norm, marker=None, color='k', lw=2, alpha=1.)
axes[0].set_xlim(0.1, 1.7)
axes[0].set_ylim(0, 30)
# axes[0].yaxis.set_visible(False)
# Second panel
xbins = np.linspace(mu_lim[0]-0.1,mu_lim[1]+0.1,contour_bins)
ybins = np.linspace(sigma_lim[0]-0.1,sigma_lim[1]+0.1,contour_bins)
axes[1].plot(sampler.flatchain[:,0], sampler.flatchain[:,1], marker=',', alpha=pt_alpha,
color='#aaaaaa', linestyle='none', zorder=-10, rasterized=True)
sigma_contour(sampler.flatchain[:,0], sampler.flatchain[:,1], bins=(xbins,ybins),
ax=axes[1], contourf=False, levels=[0.683,0.955], colors=["#000000","#555555"])
axes[1].set_xlim(*mu_lim)
axes[1].set_ylim(*sigma_lim)
axes[1].set_ylabel(r"$\sigma_{\rm WD}$ [${\rm M}_\odot$]")
# Third panel
xbins = np.linspace(-0.1,0.5,contour_bins)
ybins = np.linspace(sigma_lim[0]-0.1,sigma_lim[1]+0.1,contour_bins)
axes[2].plot(sampler.flatchain[:,2], sampler.flatchain[:,1], marker=',', alpha=pt_alpha,
color='#aaaaaa', linestyle='none', zorder=-10, rasterized=True)
sigma_contour(sampler.flatchain[:,2], sampler.flatchain[:,1], bins=(xbins,ybins),
ax=axes[2], contourf=False, levels=[0.683,0.955], colors=["#000000","#555555"])
axes[2].set_xlim(0.,0.5)
axes[2].set_ylim(*sigma_lim)
axes[2].set_ylabel(r"$\sigma_{\rm WD}$ [${\rm M}_\odot$]")
for ax in axes:
plt.setp(ax.get_xticklabels(), fontsize=14)
plt.setp(ax.get_yticklabels(), fontsize=14)
def fourth_panel(ax, data, sampler, plot_samples=False, mu_lim=(0.2,1.2), sigma_lim=(0.,1.),
pt_alpha=0.5, known_m2=None):
if known_m2 is None:
known_m2 = np.ones(len(data))*np.nan
# Fourth panel
pNSes = []
ids = np.arange(len(data['mf']))
for i in range(512):
ix = np.random.randint(len(sampler.flatchain))
pars_WD = dict(mean=sampler.flatchain[ix,0], stddev=sampler.flatchain[ix,1], bounds=(0.2,1.44))
PNS = NS_prob(pars_WD, sampler.flatchain[ix,2], data['mf'], data['m1'], known_m2)
pNSes.append(PNS[data['mf'].argsort()])
pNSes = np.array(pNSes)
bins = np.arange(len(data['mf']))
ids = np.repeat(bins[np.newaxis], 512, axis=0)
H,xedges,yedges = np.histogram2d(np.ravel(ids), np.ravel(pNSes), bins=(bins, np.linspace(0,1,50)))
ax.pcolormesh(xedges, yedges, H.T, cmap='Greys', norm=mpl.colors.LogNorm())
if np.any(data['is_NS']):
NSs = np.arange(len(data['is_NS']))[data['is_NS'][data['mf'].argsort()]]
for ns_id in NSs:
ax.plot([ns_id,ns_id], [-0.05, 0.0], linewidth=0.5, alpha=0.8,
marker=None, color='#666666')
#ax.axvline(ids[is_NS_sort].min())
buff = int(len(data['mf'])*0.05)
ax.set_xlim(-buff,len(data['mf'])+buff-1)
ax.set_ylim(-0.1,1.1)
# ax.set_ylabel(r'$P_{\rm NS}$')
In [59]:
test_data.keys()
Out[59]:
In [47]:
fig,all_axes = plt.subplots(2,4,figsize=(12,6),sharex='col')
fig2,all_axes2 = plt.subplots(2,4,figsize=(12,6),sharex='col')
for i,test_name in enumerate(test_data.keys()):
axes = all_axes[i]
data = test_data[test_name]
sampler = samplers[test_name]
pars = test_params[test_name]
n,bins,patches = axes[0].hist(data['true_m2'], color="#aaaaaa")
panel1_norm = np.sum(n * (bins[1:] - bins[:-1]))
if i == 2:
three_panel(axes, data, sampler, sigma_lim=(0.,1.5), panel1_norm=panel1_norm)
else:
three_panel(axes, data, sampler, pt_alpha=0.2, panel1_norm=panel1_norm)
axes[2].axvline(pars['f_NS'], color='k', linewidth=1., linestyle='--')
ymin,ymax = axes[0].get_ylim()
axes[0].text(0.15, ymax - 0.15*(ymax-ymin), "Test {}".format(i+1), fontsize=22, fontweight=200)
axes[0].set_ylabel("$N$")
if pars.has_key('mean_WD'):
axes[1].axhline(pars['stddev_WD'], color='k', linewidth=1., linestyle='--')
axes[1].axvline(pars['mean_WD'], color='k', linewidth=1., linestyle='--')
axes[2].axhline(pars['stddev_WD'], color='k', linewidth=1., linestyle='--')
if i == 2:
axes[0].set_xlabel(r'$M_2$ [${\rm M}_\odot$]')
axes[1].set_xlabel(r"$\mu_{\rm WD}$ [${\rm M}_\odot$]")
axes[2].set_xlabel(r'$f_{\rm NS}$')
# axes[3].set_xlabel(r'$m_f$ [${\rm M}_\odot$]')
# axes[3].set_xlabel(r'Object ID')
fourth_panel(all_axes2[i,3], data, sampler, pt_alpha=0.2)
all_axes2[i,3].set_xlabel(r'Object ID')
all_axes2[0,0].set_ylabel(r'$P_{\rm NS}$')
fig.tight_layout()
fig2.tight_layout()
fig.savefig("../documents/letter/test_12.pdf")
fig2.savefig("../documents/letter/test_3pceb.pdf")
In [60]:
slc = slice(38,43)
nm = 'GWD'
ixx = test_data[nm]['mf'].argsort()
print('sini', test_data[nm]['sini'][ixx][slc])
print('true m2', test_data[nm]['true_m2'][ixx][slc])
mff = test_data[nm]['mf'][ixx][slc].copy()
m11 = test_data[nm]['m1'][ixx][slc].copy()
m11[0] = m11[2]
mff[0] = mff[2]
print('mf', mff)
print('m1', m11)
NS_prob(dict(mean=0.724203284699, stddev=0.197681722258),
0.000740795465026,
mff, m11,
known_m2=np.array([np.nan,np.nan,np.nan,np.nan,np.nan]),
bounds_WD=(0.2,1.44))
Out[60]:
In [67]:
def four_panel(axes, data, sampler, plot_samples=False, mu_lim=(0.2,1.2), sigma_lim=(0.,1.),
pt_alpha=0.5, known_m2=np.ones(len(data))*np.nan, panel1_norm=1., contour_bins=33):
map_idx = sampler.flatlnprobability.argmax()
mean_WD,stddev_WD,f_NS = sampler.flatchain[map_idx]
print("MAP params: μ={}, σ={}, f_NS={}".format(mean_WD,stddev_WD,f_NS))
MAP_pars_WD = dict(mean=mean_WD, stddev=stddev_WD, bounds=(0.2,1.44))
MAP_f_NS = f_NS
pars_NS = dict(mean=mean_NS, stddev=stddev_NS, bounds=bounds_NS)
# First panel
m2_grid = np.linspace(0., 2., 256)
if plot_samples:
for i in range(50):
ix = np.random.randint(len(sampler.flatchain))
pars_WD = dict(mean=sampler.flatchain[ix,0], stddev=sampler.flatchain[ix,1], bounds=(0.2,1.44))
f_NS = sampler.flatchain[ix,2]
d1,d2 = model(pars_WD=pars_WD, pars_NS=pars_NS)
func = (1-f_NS)*d1.pdf(m2_grid) + f_NS*d2.pdf(m2_grid)
axes[0].plot(m2_grid, func*panel1_norm, marker=None, color='k', lw=1, alpha=0.1)
d1,d2 = model(pars_WD=MAP_pars_WD, pars_NS=pars_NS)
func = (1-MAP_f_NS)*d1.pdf(m2_grid) + MAP_f_NS*d2.pdf(m2_grid)
axes[0].plot(m2_grid, func*panel1_norm, marker=None, color='k', lw=2, alpha=1.)
axes[0].set_xlim(0.1, 1.7)
axes[0].set_ylabel("$N$")
# Second panel
xbins = np.linspace(mu_lim[0]-0.1,mu_lim[1]+0.1,contour_bins)
ybins = np.linspace(sigma_lim[0]-0.1,sigma_lim[1]+0.1,contour_bins)
axes[1].plot(sampler.flatchain[:,0], sampler.flatchain[:,1], marker=',', alpha=pt_alpha,
color='#aaaaaa', linestyle='none', zorder=-10, rasterized=True)
sigma_contour(sampler.flatchain[:,0], sampler.flatchain[:,1], bins=(xbins,ybins),
ax=axes[1], contourf=False, levels=[0.683,0.955], colors=["#000000","#555555"])
axes[1].set_xlim(*mu_lim)
axes[1].set_ylim(*sigma_lim)
axes[1].set_ylabel(r"$\sigma_{\rm WD}$ [${\rm M}_\odot$]")
# Third panel
xbins = np.linspace(-0.1,0.5,contour_bins)
ybins = np.linspace(sigma_lim[0]-0.1,sigma_lim[1]+0.1,contour_bins)
axes[2].plot(sampler.flatchain[:,2], sampler.flatchain[:,1], marker=',', alpha=pt_alpha,
color='#aaaaaa', linestyle='none', zorder=-10, rasterized=True)
sigma_contour(sampler.flatchain[:,2], sampler.flatchain[:,1], bins=(xbins,ybins),
ax=axes[2], contourf=False, levels=[0.683,0.955], colors=["#000000","#555555"])
axes[2].set_xlim(0.,0.5)
axes[2].set_ylim(*sigma_lim)
axes[2].set_ylabel(r"$\sigma_{\rm WD}$ [${\rm M}_\odot$]")
# Fourth panel
pNSes = []
ids = np.arange(len(data['mf']))
for i in range(512):
ix = np.random.randint(len(sampler.flatchain))
pars_WD = dict(mean=sampler.flatchain[ix,0], stddev=sampler.flatchain[ix,1], bounds=(0.2,1.44))
PNS = NS_prob(pars_WD, sampler.flatchain[ix,2], data['mf'], data['m1'], known_m2)
pNSes.append(PNS[data['mf'].argsort()])
pNSes = np.array(pNSes)
bins = np.arange(len(data['mf']))
ids = np.repeat(bins[np.newaxis], 512, axis=0)
H,xedges,yedges = np.histogram2d(np.ravel(ids), np.ravel(pNSes), bins=(bins, np.linspace(0,1,50)))
axes[3].pcolormesh(xedges, yedges, H.T, cmap='Greys', norm=mpl.colors.LogNorm())
try:
if np.any(data['is_NS']):
NSs = np.arange(len(data['is_NS']))[data['is_NS'][data['mf'].argsort()]]
for ns_id in NSs:
axes[3].plot([ns_id,ns_id], [-0.05, 0.0], linewidth=1., alpha=0.8,
marker=None, color='#333333')
except:
pass
buff = int(len(data['mf'])*0.05)
axes[3].set_xlim(-buff,len(data['mf'])+buff-1)
axes[3].set_ylim(-0.1,1.1)
axes[3].set_ylabel(r'$P_{\rm NS}$')
for ax in axes:
plt.setp(ax.get_xticklabels(), fontsize=14)
plt.setp(ax.get_yticklabels(), fontsize=14)
In [63]:
fig,axes = plt.subplots(2,4,figsize=(16,9))
# Test 1 and 2
for i,test_name in enumerate(['GWD','GWD+GNS']):
pars = test_params[test_name]
data = test_data[test_name]
row = axes[i]
n,bins,patches = row[0].hist(data['true_m2'], color="#aaaaaa")
panel1_norm = np.sum(n * (bins[1:] - bins[:-1]))
four_panel(row, data, samplers[test_name], plot_samples=False,
sigma_lim=(0.,1.5), panel1_norm=panel1_norm, contour_bins=39)
row[2].axvline(pars['f_NS'], color='k', linewidth=1., linestyle='--')
ymin,ymax = row[0].get_ylim()
row[0].text(0.15, ymax - 0.15*(ymax-ymin), "Test {}".format(i+1), fontsize=22, fontweight=200)
row[1].axhline(pars['stddev_WD'], color='k', linewidth=1., linestyle='--')
row[1].axvline(pars['mean_WD'], color='k', linewidth=1., linestyle='--')
row[2].axhline(pars['stddev_WD'], color='k', linewidth=1., linestyle='--')
row[0].set_xlabel(r'$M_2$ [${\rm M}_\odot$]')
row[1].set_xlabel(r"$\mu_{\rm WD}$ [${\rm M}_\odot$]")
row[2].set_xlabel(r'$f_{\rm NS}$')
row[3].set_xlabel(r'Object ID')
fig.tight_layout()
fig.savefig("../documents/letter/test_1_2.pdf")
In [68]:
fig,axes = plt.subplots(2,4,figsize=(16,9))
# Test 3
test_name = 'UWD+GNS'
pars = test_params[test_name]
data = test_data[test_name]
row = axes[0]
n,bins,patches = row[0].hist(data['true_m2'], color="#aaaaaa")
panel1_norm = np.sum(n * (bins[1:] - bins[:-1]))
four_panel(row, data, samplers[test_name], plot_samples=False,
sigma_lim=(0.,1.5), panel1_norm=panel1_norm, contour_bins=31)
row[2].axvline(pars['f_NS'], color='k', linewidth=1., linestyle='--')
ymin,ymax = row[0].get_ylim()
row[0].text(0.15, ymax - 0.15*(ymax-ymin), "Test 3", fontsize=22, fontweight=200)
# PCEB
row = axes[1]
n,bins,patches = row[0].hist(pceb_MD['M_wd'][~np.isnan(pceb_MD['M_wd'])], color='#bbbbbb')
ymin,ymax = row[0].get_ylim()
row[0].text(0.15, ymax - 0.15*(ymax-ymin), "PCEB".format(i+1), fontsize=22, fontweight=200)
panel1_norm = np.sum(n * (bins[1:] - bins[:-1]))
four_panel(row, pceb_MD, sampler_real_data_pceb, plot_samples=False,
panel1_norm=panel1_norm, contour_bins=39)
row[0].set_xlabel(r'$M_2$ [${\rm M}_\odot$]')
row[1].set_xlabel(r"$\mu_{\rm WD}$ [${\rm M}_\odot$]")
row[2].set_xlabel(r'$f_{\rm NS}$')
row[3].set_xlabel(r'Object ID')
fig.tight_layout()
fig.savefig("../documents/letter/test_3_pceb.pdf")
In [69]:
fig,axes = plt.subplots(1,4,figsize=(16,4.5))
four_panel(axes, elm_WD, sampler_real_data, plot_samples=True, known_m2=known_m2)
ymin,ymax = axes[0].get_ylim()
axes[0].text(0.15, ymax - 0.15*(ymax-ymin), "ELM".format(i+1), fontsize=22, fontweight=200)
axes[0].set_ylabel(r"$p(M_2)$")
axes[0].set_xlabel(r'$M_2$ [${\rm M}_\odot$]')
axes[1].set_xlabel(r"$\mu_{\rm WD}$ [${\rm M}_\odot$]")
axes[2].set_xlabel(r'$f_{\rm NS}$')
axes[3].set_xlabel(r'Object ID')
fig.tight_layout()
fig.savefig("../documents/letter/real-data.pdf")
In [872]:
np.savetxt("posterior_samples.txt", sampler_real_data.flatchain, delimiter=",", header="mu_WD, sigma_WD, f_NS")
In [ ]: